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TECHNICAL MEMORANDUM 


I. INTRODUCTION 


The purpose of this report is to present the NASA-Marshall Space Flight Center 
solar activity long-range statistical estimation technique and the basis for the proce- 
dure used. The motivation for issuing the information monthly is the need for a 
current and systematic input for use in an upper atmosphere der '•^ty model employed 
in satellite orbital lifetime estimates. The mission analysis and planning for spacecraft 
operations require estimates of orbital lifetimes for launches up to 15 years in the 
future relative to various orbital altitudes, inclinations, and eccentricities. 


The upper atmosphere density model used [1] is based on the work an models 
developed by Jacchia [2,3] and his colleagues at the Smithsonian Astrophysical Obser- 
vatory (SAO). A more recent model [4] uses essentially the same type solar activity 
input and gives more detailed information on number densities for the various constit- 
uents. To compute neutral density at the pertinent orbital altitudes and times, the 
model requires as inputs the 162-day mean and daily value of 10.7 cm solar flux 7) 


and the 3-hour interval geomagpietic index (A^) 6 to 9 hours prior to the time of the 


desired density calculation. The variability of the solar flux and geomagrietic activity 
precludes the long-term prediction of these daily and 3-hourly values with ary degree 
of accuracy; therefore, smoothed values of these indices, estimated using a least- 
square linear regression procedure, are employed as inputs to the model. The model 
then predicts the neutral density values required in the computation of satellite orbital 
lifetime estimates. 


This study presents the results of comparison studies conducted in the Atmos- 
pheric Sciences Division, Space Sciences Laboratory, Marshall Space Flight Center 
(MSFC) , on various ap plica dons of statistical prediction methods for short-term 
(months) and long-term (years) forecasting of solar activity. This is one of many 
approaches to the estimation of future solar flux and geomagnetic activity prepared by 
those who have and are now doing research on the subject. There is no generally 
accepted statistical or deterministic procedure based on any known physical mechanism 
by the scientific community. The one reported here is an attempt to improve upon the 
statistical confidence of a linear digression procedure developed around 1949 by McNish 
and Lincoln. 


The basic theory first suggested and applied by McNish and Lincoln [5] and 
later modified by Boykin and Ricliards [6j and improved by Rowe and Avaritt f^] is 
gi\'en in Appendix A. An excellent survey is also presented by Scissum [8] and 
Slutz, et al. [9]. 


II. TECHNIQUE 


A Historical Background 

McNish and Lincoln [5] suggested that the estimation of a future solar cycle, 
based on the mean approximation of all past solar cycles, could be improved by adding 
to the mean a correction proportional to the departure of the current values of the 
cycle from the mean cycle. Howcsver, this technique was not recommended for making 
projections longer than 1 year in the future. 



f ir\ 
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Yu. I. Vitinskii [10] conducted a very extensive survey and analysis of solar 
activity forecasting methods. While recognizing the considerable importance of the 
problem and encouraging studies of active processes taking place on the Sun as a con- 
tributiixi toward the solution of the proolem , he well stated what is still the current 
situation. Vitinskii concluded: *'...we have shown that the reliability of the results 
obtained using these methods still leaves much to be desired." However, his analysis 
showed that the linesu* regression method gives the best acciiracy for future solar 
activity estimates up to a year in advance. For estimates several years in advance, 
the linear regression method becomes increasingly less accurate. The uncertainty of 
tlie linear regression method is given by the standard deviation estimat;:s as computed 
by the least-square linear regression method reliability procedure. With the lack of 
any proven deterministic prediction scheme based on physical laws, the statistical 
predictions produced by this technique are used as a basis for calculating upper 
atmospheric density values needed for mission planning and analysis. We thereby 
have a best estimate and confidence bounds bat x on the available data record of 
solar activity. 

Using data from two additional solar cycles, Boykin [6] modified the McNish- 
Lincoln technique so that the Zurich- smoothed sunspot number (R) could be estimated 
for 10 years in advance, and at quai*terly rather than yearly intervals. 

Since the measured Fjq y data base is relatively short (1947 to date), the data 

base was extended back to 1749 »*sing the following equation^ to convert recorded 
smoothed (R) data to smoothed solar flux data: 

Ffo y = 49. 4 + 0.97 R + 17.6 exp (-0.035 R) . (1) 

The recorded data for Fjq y and R from 1947 to the present were used to 

derive equation (1). The linear correlation coeffi<dent_for eqmtion (1) is 0.98. 
Equation (1) is used to construct smoothed solar flux Fjq y data for the period 1749 

to 1946. After 1946, the measured values of the daily solar flux F^q ^ are used to 

compute the ^0. This combined data base was then stored by cycle number and 

point number, Lagrange interpolated and inputed into the main linear regression 
subroutine. 

All solar flux cycles used in the regression model data base prior to 1947 use 
the Zurich established sunspot number maxima as the starting point. After 1947, the 
maxima are based on the actual maxima from the smoothed solar flux values computed 
from the measured F^q ^ values. 

The output of the solar activity long-range statistical scheme is an estimate of 
the 97.7, 50, and 2,3 percentile values of the F.j. „ arranged in a monthly time 
sequence. 


1. Developed in collaboration with Jack Slowey of Smithsonian Astrophysical 
Observatory. 


Since the measured geomagnetic data base relatively short (1932 to date), 

the data base was extended back to 1884 using equation (2) to convert geomagnetic 

C data to geomagnetic A data. 

* P 

Ap = Exp (0.8267 + 4.3661 C. - 4.0832 C.^ + 2.1795 C.® - 0.3536 C.^ . (2) 

The recorded data for Ap and C. from 1932 to the present were used to derive 

equation (2). The linear correlation coefficient for eq\iation (2) is 0.99. Cycles 13 
through 16 computed Ap valu'^s were then combined with actual Ap values from 1932 

through present. This combined data base was then smoothed and stored by cycle 
number and point number, Lagrange interpolated and put into the main linear regres- 
sion subroutine. 

The output of the solar activity long-range statistical scher'3 is an estimate of 
the 97.7, 50, and 2.3 percentile values of Ap arranged in a mont'nly time sequence. 


6. Lagrangian Least Squares Method 

The Lagrangian least-squares prediction method outlined in the following steps 
was compared with the conventional application of the McNish- Lincoln method and 
also with various applications of the Lagrangian method to determine the methods 
which gave the best predictions of past cycles of solar activity. 

1) Select all maximum (or minimum) points. 

2) Determine all the Periods {P.}. 

3) Determine the average period P. 

4) Divide all cycle periods (in step 2) by P (the average period) . 

This will give an average time increment Tj s P^/P to be used for the jth cycle. Note 
that Tj will be less thsin, equal to, or greater than 1 mor Lh. 


5) Consider each cycle a block of data by selecting m = P sunspot (or flux) 
points in each .-ycle. This will require some interpolation to obtain each sunspot (or 
flux) corresponding co all m points per cycle, since points are not being selected at 
the 1-month increment. 


6) When all the cycles have been blocked in this manner, the McNish- Lincoln 
metliod can be applied to produce a prediction of the mean and ±lo, 2a, and 3o values 
of Fio rj for the next cycle, assuming that the distribution of values is Gaussian. 

7) The predicted cycle time increment is approximated by averaging all the past 
cycle time increments. It just so happens that this time increment will be exactly 
equal to 1. 
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(3) 


N is the number of past cycles. 

P. is the period of the jth cycles. 

8) The interpolation mentioned in step 5 is accomplished by applying a conven- 
tional Lagrangian method from the MSFC math-pack computer program library [11]. 


9) The predicted cycle period will be the average period P plus or minus its 

2o ulue; i.e.. P ± 2o . 

P 

10) The 2o error band on the time axis combined with the 2o_ error band on 
p R 

the sunspot (or flux) values will give the "error space" on the two-dimensional plot 

of R versus t. 


III. SINGLE CYCLE PREDICTION 


The preceding procedure was applied, with the number of points per cycle (m) 
being 132. This number was selected because it is the closest integer to the average 
period in months. A single cycle was then predicted using the previous 16 cycles 
(from 1749) in order to predict cycle 17 for comparison with the actual observations. 
This process was then repeated in order to predict cycles 17, 18, 19, and 20 by the 
conventions’ McNish- Lincoln method and the Lagrangian least-squares approach. 


This particular statistic gives a relatively simple and easily applied criterion for 
judging and comparing goodness of fit of various models (or expected values) of 
observed uata. I’or large degrees of freedom, the standard chi square disiribution 
approaches the normal distribution. It also gives a standard, acceptable and fair way 
of determining those r..odels (or portions of the models) wliich fall within the con 
ventional 95 percentile level of confidence. 


Table 1 gives a summary of the results of the preceding ..iethod. The entries 

O 

are the calculated (0.05) by Pearson's chi-square statistics [12]; i.e.. 


2 

X = 


m 


(O, 


¥ 


where O. is observed and E. is the expected (predicted) value. 


( 4 ) 
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TABLE 1. SINGLE CYCLE x^(0.05) 


Method 

17 

18 

19 

20 

McNish- Lincoln 





1 Coefficient 

721 

1687 

268 

400 

4 Coefficient 

1404 

6128 

236 

304 

Lagrangian 





1 Coefficient 

277 

1193 

155 

177 

4 Coefficient 

753 

513 

1221 

183 


2 2 

Note: These are to be compared with the \ (0.05, 127) and x (0.05, 130) 

values from the standard chi-square tables, which are 154 and 158. 


Table 1 shows that the Lag^rangian least squares approach produces a signifi- 
cantly better prediction of the next solar cycle than the McNish- Lincoln technique 
with the exception of cycle 18. 


IV. EXTENDED AND DOUBLE CYCLE PREDICTION 


This section presents descriptions of the six methods for extended predictions 
which were compared, again using Pearson's chi-square statistics [see Eq. (4)]. 

The double-cycle McNish -Lincoln technique applied the standard regression 
program to the data blocked in 262 points per block, which is approximately 22 years 
or two average cycle periods. The 262 paints started at the maximum and ended 262 
rix>nths later. 

The double-cycle Lagrangian technique blocked the data into double cycles of 262 
points which started at the maximum and ended at a maximum two cycles later. This 
necessitates interpolation because the time increment between points varies from block 
to block. 

The technique of Langfrangiun prediction connected to the mean at the minimum 
blocks the data into single cycles of 132 points per cycle. The 132 points start at the 
maximum and end at the next maximum. A future cycle is predicted to the next minimum 
and then patched (connected with a cubic at the nearest inflection point 2 years before 
the maximum) to the mean cycle for the remaining part of the extended prediction (see 
Appendix B for cubic connection). 

The technique of Lagrangian prediction connected to the mean at the maximum 
technique is similar to the preceding technique except that the prediction is over a full 
future cycle to the next maximum and then connected to the mean cycle for the remain- 
ing part of the extended prediction. 
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The iterated Lagrangian technique predicts over a full cycle, treats the predicted 
data as actual data, and then predicts over another future cycle. 

The Lagrangian derivative technique includes the influence of the rate of change 
of the deviations from the mean cycle in the regression scheme. The regression model 
has this term added in the same way as the deviations are included. 


m 

m+k,j m+k,j . q,k q.j q,k q, 




(5) 


If the derivative influence is zero, the d . coefficients will be zero or the first m Aft's 
will be zero. 


Each of the preceding methods was applied to predict cycles 15 through 20 by 
using all the previous cycles to the one being predicted. Since actual or observed 
data and expected or predicted data were then available. Pearson's [12] chi-square 
could be calculated at each point of the cycle (beyond, say, 10 or 20 for a significant 
number of degrees of freedom). The results of applying the methods above are shown 
in Figures I through 40. 

The even figure plots (2, 4, 6, .... 40) give the predicted solar flux (.), the 
actual smoothed values (♦), the plus 2 o error (+) , and the minus 2o (-) error bands. 
The plus and minus 2o error was calculated using the predicted (expected) values of 
the current cycle or cycles minus the actual 13-month smoothed Lagrangian values at 
corresponding points of all previous cycles for the deviations. 


1 



where i refers to the Lagrangian point in the cycle and N is the number of known past 
cycles. 

2 

The odd figure plots (1. 3, 5 .... 39) give the \ versus point number (i) in the 

*> 

cycle together with the standard \" (0.05) at the 0.05 level of significance versus 
points (i) in the cycle. Table 2 presents tho overall chi-square comparisons for the 
full extended predictions. It can be seen again that the Lagranguui technique gave 
a better chi-square than the McNish- Lincoln, but it still was not within tho 95 per- 
centile level of confidence given in the right column, except for cycles 15 and 16. 

It may also be noted that not one of the metlwds tried did well on cycles 18 and 19. 
Tlus suggests that the long period effect (7 or 8 cycles) referred to in the literature 
and in Section VI cannot be predicted by any of the metlxids tried thus far. From 
the plots one can d 'termine 1k>w well tho prediction metliod is doing as it progresses 
through the cycle and beyond. The plots also sliow when the metlvid is within tl>e 
0.05 level of significance. 




Figure 1. Chi square versus point number (nominally months), 13-month 
r.agrangian-smoothed flux (extended by median cycle from min.), 1 
coefficient linear model using 1 through 14 predicting cycles 15, 

and 16 from point 1. 



Fil virc 2. A^redicted and actual flux versus time, 13-month La gran gian- smoothed 
flux (extended oy median cycle from min.), 1 coefficient linear model using 
cycles 1 through 14 predicting cycles, 15 and 16 from point 1. 
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Figure 3. Chi square versus point number (nominally months) 13-month 
La gran gian- smoothed flux (extended by median cycle from min.) 1 
coefficient linear model using cycles 1 through 15 predicting cycles 

16 and 17 from point 1. 



Figure 4. Predicted and actual flux versus time 13-month Lagrangian- smoothed 
flux (extended by median cycle from min.) 1 coefficient linear model using 
cycles i through 15 predicting cycles 16 and 17 from point 1. 






Figure 5. Chi square versus point number (nominally months) 13-month 
Lagrangian- smoothed flux (extended by median cycle from min.) 1 
coefficient linear model using cycles 1 through 16 predicting cycles 

17 and 18 from point 1. 



Figure 6. Predicted and actual flux versus time 13-month Lagrangian smoothed 
flux (extended by median cycle from min.) 1 coefficient linear model using 
cycles 1 through 16 predicting 17 and 18 from point 1. 




Figure 7. Chi square versus point number (nominally months), 13-month Lagrangian- 
smoothed flux (extended by median cycle from min.) 1 coefficient linear model using 
cycles 1 through 17 predicting cycles 18 and 19 from point 1. 



Figure 8. Predicted and actual flux versus time 13-month Lagrangian- smoothed 
flux (extended by median cycle from min.) 1 coefficient linear model using 
cycles 1 through 17 predicting cycles 18 and 19 from point 1. 
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Figure 9. Chi square versus point number (nominally months) 13-month 
Lagrangian- smoothed flux (extendea median cycle from min.) 1 
coefficient linesu* model using cycles 1 through 18 predicting cycles 

19 and 20 from point 1. 
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Figure 10. Predicted and actual flux versus time 13-month La gran gfian- smoothed 
flux (extended by median cycle from min.) 1 coefficient linear model using 
cycles 1 through 18 predicting cycles 19 and 20 fl*om point 1. 
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Figure 13. Chi square versus point number (nominally months) 13-month 
Lagrangian- smoothed flux (extended by median cycle from max.) 1 
coefficient linear model using cycles 1 through 15 predicting cycles 

16 and 17 from point 1. 



Figure 14. Predicted and actual flux versus time 13-month Lagrangian-smoothed 
flux (extended by median cycle from max.) 1 coefficient linear model using 
cycles 1 through 15 predicting cycles 16 and 17 from point 1. 
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Figure 15. Chi square versus point number (nominally months) 13-month 
Lagrangian-smoothed flux (extended by median cycle from max.) 1 
coefficient linear model using cycles 1 through 16 predicting cycles 

17 and 18 from point 1. 



VtAII 


Figure 16. Predicted and actual flux versus time 13- month Lagrangian-smoothed 
flux (extended by median cycle from max.) 1 coefficient linear model using 
cycles 1 through 16 predicting cycles 17 and 18 from point 1. 



Figure 17. Chi square versus point number (nominally months) 13-month 
La gran gian- smoothed flux (extended by median cycle from max.) 1 
coefficient linear model using cycles 1 through 17 predicting cycles 

18 and 19 from point 1. 



YEAR 

Figure 18. Predicted and actual flux versus time 13-month Lagrangian- smoothed 
flux (extended by median cycle from max.) 1 coefficient linear model using 
cycles 1 through 17 predicting cycles 18 and 19 from point 1. 
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Figure 19. Chi square versus point number 13-month Lagrangian -smoothed flux 
(extended by median cycle ‘'"cm max.) 1 coefficient linear model using 
cycles 1 through 18 predicting cycles 19 and 20 from point 1. 



Figure 20. Predicted and actual flux versus time 13-month La gran gian- smoothed 
flux (extended by median cycle from max.) 1 coefficient linear model using 
cycles 1 through 18 predicting cycles 19 and 20 from point 1. 
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Figure 21. Chi square versus point number (noninally months) 13-month 
La gran gian- smoothed flux (extended by iterated McNish-Lincoln) 1 
coefficient linear model using cycles 1 through 14 predicting cycles 

15 and 16 from point 1. 



Figure 22. Predicted and actual flux versus time 13-month Lagrangian- 
smoothed flux (extended by iterated McNish-Linooln) 1 
coefficient linear model using cycles 1 through 14 predicting 
cycles 15 and 16 from point 1. 



ORIGINAL PAGE IS 
OF POOR QUALITY 



Figure 23. Chi square versus point number (nominally months) 13-month 
La gran gian- smoothed flux (extended by iterated McNish-Lincoln) 1 
coefficient linear model using cycles 1 through 15 predicting 
cycles 16 and 17 from point 1. 



Figure 24. Predicted and actual flux versus time 13-month La gran gian -smoothed 
flux (extended by iterated McNish-Linooln) 1 coefficient linear model using 
cycles 1 through 15 predicting cycles 16 and 17 from point 1. 
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Figure 25. Chi square versus point number (nominally months) 13- month 
Lagrangian-sn)oothed flux (extended by iterated McNish-Lincoln) 1 
coefficient linear model using cycles 1 through 16 predicting 
cycles 17 and 18 from point 1. 



Figure 26. Predicted and actual flux versus time 13-month Lagrangian- smoothed 
flux (extended by iterated McNish-Lincoln) 1 coefficient linear model using 
cycles 1 through 16 predicting cycles 17 and 18 from point 1. 
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Figure 27. Chi square versus point number (nominally months) 13-month 
Lagrangian -smoothed flux (extended by iterated McNish-Lincoln) 1 
coefficient linear model using cycles 1 through 17 predicting 
cycles 18 and 19 from point 1. 



VIAR 

Figure 28. Predicted and actual flux versus time 13 month Lagrangian-smoothed 
flux (extended by iterated McNish-Lincoln) 1 coefficient linear model using 
cycles 1 through 17 predicting cycles 18 and 19 from point 1. 
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Figure 29. Chi square versus point number (nominally months) 13-month 
La gran gian- smoothed flux (extended by iterated McNish-Lincoln) 1 
coefficient linear model using cycles 1 through 18 predicting 
cycles 19 and 20 from point 1. 
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Fig^e 30. Predicted and actual flux versus time 13-month Lagrangian-smoothed 
flux (extended by iterated McNish-Linooln) 1 coefficient linear model using 
cycles 1 through 18 predicting cycles 19 and 20 from point 1. 
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Figure 31. Chi square versus poL t number (nominally months) 13-month 
) Lagrangian- smoothed flux (extended by derivative and connected to 
median cycle from max. ) 1 coefficient linear model using cycles 1 
through 14 predicting cycles 15 and 16 from point 2. 



Figure 32. Predicted and actual flux versus time 13-month Lagrangian -smoothed 
flux (extended by derivative and connected to median cycle from max.) 

1 coefficient linear model using cycles 1 through 14 predicting 
cycles 15 and 16 from point 2. 
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Figure 33. Cni square versus point number (nominally months) 13- month 
Lagrangian-smoothed flux (extended by derivative and connected to median 
cycle from max.) 1 coefficient linear model using cycles 1 through 15 
predicting cycles 16 and 17 from point 2. 
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Figfure 34. Predicted and actual flux versus time 13-month Lagrangian-smoothed 
flux (extended by derivative and connected to median cycle from max.) 1 
coefficient linear model using cycles 1 through 15 
predicting cycles 16 and 17 from point 2. 
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Figure 35. Chi t'quare versus point number (nominally months) 13-month 
Lagrangian-smoothrid fl\ix (extended by derivative and connected to median 
cycle from max.) 1 coefficiert linear model using cycles I through 16 
predicting cycles 17 and 18 from point 2. 



Figure 36. Predicted and actual flux versus time 13-month Lagrangian-smoothed 
flux (extended by derivative and connected to median cycle from max.) 1 
coefficient linear model using cycles 1 through 16 
predicting cycles 17 and 18 from point 2. 
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Figure 37. Chi square versus point number (nominally months) 13-month 
Lagrangian-srooothed flux (extended by derivative and connected to mecUan 
cycle from max. ) 1 coefficient linear model using cycles 1 through 17 
predicting cycles 18 and 19 from point 2. 



Figure 38. Predicted and actual flux versus dme 13-month La gran gian- smoothed 
flux (extended by derivative and connected to median cycle from max.) 1 
coefficient linear model using cycles 1 through 17 
predicting cycles 18 and 19 from point 2. 
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Figure 39. Chi square versus point number (nominally months) 13-month 
Lagrangian- smoothed flvix (extended by derivative and connected to median 
cycle firom max.) 1 coefficient linear model using cycles 1 through 18 
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Figure 40. Predicted and actual flux versus time 13-month Lagrangian-smoothed 
flux (extended by derivative and connected to median cycle from max.) 1 
coefficient linear model using cycles 1 through 18 
predicting cycles 19 and 20 from point 2. 
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TABLE 2. EXTENDED CYCLE 


Method 

15 ft 16 

16 ft 17 

17 ft 18 

18 ft 19 

19 ft 20 

X^(0.05) 

Double -cycle 







1. McNish-Lincoln 

1891 


7597 


738 

300 

From Point 12 

1828 


7753 


1000 

287 

From Point 24 



8179 


912 

274 

Double-cycle 







2. Lagrangiw 

190 


4000 


650 

300 

From point 12 

204 


3521 


643 

287 

From point 24 

213 


3075 


383 

274 

Lagrangian 







3. Patched at min. 

214 

426 

1427 

3306 

454 

242 

From point 12 

175 

425 

1344 

3371 

514 

229 

From point 24 

181 

439 

1051 

3250 

459 

216 

Lagrangian 







4. Patched at max. 

238 

351 

1348 

3074 

479 

242 

From point 12 

218 

355 

1243 

3104 

597 

229 

From point 24 

227 

362 

1155 

2549 

634 

216 

Lagrangian 







5. Iterated 

248 

676 

2562 

2406 

348 

242 

From point 12 

207 

728 

2216 

2355 

686 

229 

From point 24 

231 

653 

1589 

725 

1038 

216 

Lagrangian 







6. Derivative 

244 

381 

1430 

3057 

617 

242 

From point 12 

255 

344 

1312 

2804 

727 

229 

FrcMD point 24 

248 

340 

1249 

2714 

435 

216 


V. CONSTRUCTION OF MODEL DATA BASE 


The NASA -Marshall Space Flight Center solar activity long-range statistical estima- 
tion technique utilizes the 13-month Zurich procedure smoothed values of solar flux 
(^10 7) the data base. For the period from August 1947 to the present, the actual 

measured values of solar flux, smoothed by the Ziirich procedure, are used in the pro- 
gram. Prior to 1947, the solar flux was estimated using available Zurich sunspot numbers 
according to the procedure presented in Section II of this report. The estimated data 
plus actual measured data are then used in the main linear regression program to predict 
the remainder of the current cycle (peak-to-peak or minimum -to-minimum depending on 
initialization point), including the statistical confldence bounds, 97.7 and 2.3 percentile 
values. Table 3 provides a summary of the data base characteristics. 
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TABLE 3. FLUX 


Cycle 

Year 

of 

Minimum 

Year 

of 

Maximum 

Rise 
Time 
In Years 

Fall 
Time 
In Years 

Min 

to 

Min 

Period 

Ma:c 

to 

Max 

Period 

0 

1745 

1750.3 

5.3 

4.9 

10.2 

11.2 

1 

1755.2 

1761.5 

6.3 

5.0 

11.3 

8.2 

2 

1766.5 

1769.7 

3.2 

5.8 

9.0 

8.7 

3 

1775.5 

1778.4 

2.9 

6.3 

9.2 

9.7 

4 

1784.7 

1783.1 

3.4 

10.2 

13.6 

17.1 

5 

1798.3 

18C5.2 

6.9 

5.4 

12.3 

11.2 

6 

1810.6 

1816.4 

5.8 

6.9 

12.7 

13.5 

7 

1823.3 

1829.9 

6.6 

4.0 

10.6 

7.3 

8 

1833.9 

1837.2 

3.3 

6.3 

9.6 

10.9 

9 

1843.5 

1848.1 

4.6 

7.9 

12.5 

12.0 

10 

1856.0 

1860.1 

4.1 

7.1 

11.2 

10.5 

11 

1867.2 

1870.6 

3.4 

8.3 

11.7 

13.3 

12 

1878.9 

1883.9 

5.0 

5.7 

10.7 

10.2 

13 

1889.6 

1894.1 

4.5 

7.6 

12.1 

12.9 

14 

1901.7 

1907.0 

5.3 

6.6 

11,9 

10.6 

15 

1913.6 

1917.6 

4.0 

6.0 

10.0 

10.8 

16 

1923.6 

1928.4 

4.8 

5.4 

10.2 

9.0 

17 

1933.8 

1937.4 

3.6 

6.7 

10.3 

10.3 

18 

1944.1 

1947.7 

3.6 

6.5 

10.1 

10.5 

19 

1954.2 

1958.2 

4.0 

6.4 

10.4 

12.4 

20 

21 

1964.6 

1976.4 

1970.6 

6.0 

5.8 

11.8 



STD . Deviation a 

1.2 

1.34 

1.23 

2.2 

21 Cycle Average 
20 Cycle Average 

4.6 

6.4 

11.02 

11.02 

Interpolation Points 

53 

79 




VI. SMOOTHED PERIOD AND AMPLITUDE VARIATION 


Using the sunspot periods obtained from Waldmeier [13] , Chernosky and Hagan 
[14] , or Schove [15] , Table 3 was constructed up to 1947. Beyond 1947, actual flux 
observations were used. A five-cycle smoothing of these periods, rise times and fall 
times versus cycle number, is given in Figures 41 and 42. Figure 43 gives the five- 
cycle smoothed maximum amplitudes. These curves suggest a seven- to eight-cycle 
variation in the periods as well as the maximum ampliti.des. 


To further study this variation , autocorrelation and power spectrum analyses 
were run. Figures 44 and 45 give the results. It should be pointed out that these 
were determined from the Lagrangian interpolated data. This explains the rather 
regular variation at 132-point intervals on the autocorrelation and the sharp spike at 
1/132 frequency on the power spectrum. 
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Figure 41. 5*cycle smooth periods (min. to nrdn. ) , 0 cycle = year 1745 

min and year 1750.3 max. 
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Figure 43. 5 -cycle smoothed sunspot atrplitudti of IS^month 
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Figure 45. Power spectrum of cycles 1 through 20 Lagrangian-sroooth flux. 

Even though the Lagrangian interpolated data forced the preceding variations, it 
did not force the long-term variation (eight cycles) in the autocorrelation or the lower 
frequency peak at 0.0011 in the power spectrum. This gives some evidence to support 
the application of a connection to the average period for future predictions. 

The following examples use the apparent systematic variation in the periods in 
Figures 41 and 42 to project the maximum to minimum fall time for cycle 21 and the 
maximum to maximum period for the same cycle . 


Maximum to Minimum Fall Time— Cycle 21: 


From Figure 41 the five -cycle smoothed fall time for cycle 19 is - 6. 5 years. Now, 
if the last four fall times for cycles 17, 18, 19, and 20 plus the unknown fall time for 
cycle 21 are averaged and equated to the extrapolated smoothed fall time (7 years) . 
we can then solve for the unknown fall time for cycle 21. 


1 

5 


(F 


17 


+ F 


18 


"^ 19 "^ 


20 





(7) 
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where 


= 6.7 years (maximum to minimum period for cycle 17) 

Fig = 6.5 years 
Fjg = 6. 4 years 

^20 “ ® y®®^® 

F 21 = Unknown 

^19 “ ^ years (five-cycle smoothed period at midpoint of smoothing interval). 
From Eq. (7) one can now solve for Fg^; i.e. . 

^21 ~ ^^19 ‘ ^^17 *^18 *^19 * *^ 20 ^ 


= 5(6. 5) - (6. 7 + 6. 5 + 6. 4 + 5. 8) , 


or 


^21 ~ years maximum to minimum. 


Minimum to Maximum Rise Time -Cycle 21-22: 

Applying the same technique to the minimum to maximum rise times, the follow- 
ing results are obtained: 


R 


17 


3.6 


Ri8 = 3.6 


R 


19 


4.0 


R 


20 


6.0 


R 2 J = Unknown 
^9= 4.5 


*^21 ®^^ 9 ^ ' ^*^17 ^ *^18 * *^19 * ** 20 ^ 

= 5(4. t»; - (3.6 + 3.6 + 4.0 + 6.0), 

R 2 J = 5.3 years minimum to maximum period. 
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VII. COMPUTER PROGRAMS 


The computer programs to calculate the prv ;ted cycles and perform sub- 
sequent statistical analysis were developed by Co«nputer Sciences Corporation under 
contract NAS 8-31640. The programs were implemented in FORTRAN on the UNIVAC 
1100/80 computer at Marshall Space Flig' t Centt^r. 

Although six distinct program versions were used to obtain the numerical 
results of Table 2, all programs wore essentially variations on the single theme of 
linear regression. Gaussian reduction was used to solve the normal equations via 
subroutine GJR , contained in the UNIVAC "Math-Pack" subroutine library. The 

2 

evaluation of the inverse \ function was accomplished through subroutine CHIN 
(also drawn from the "math -pack") and through normal approximation for large 
degrees -of- freedom. Graphic output was performed off-line on an FR-80 micro- 
graphics unit. 

Space prohibits including listings for all six programs. However. FORTRAN 
source code for the liagrangiiui prediction connected at cycle maximum appears in 
Appeiuiix C. All of the system subroutines are not listed but, rather, only the main 
program and the McNish- Lincoln subroutine described in Appendix A. 

VIll. APPLICATIONS FOR SPACE PROGRAM DESIGN 
AND MISSION ANALYSIS 

The solar flux and geomagnetic activity arc principal inputs to the upvKir atmos- 
phere models and thus the ambient density values computed for use in spacecraft 
orbital lifetime estimation contml analyses programs. For those spacecraft projects 
wliich rotiuii*e a design lifetime at a given orbital altitude(s) and/or a specified con- 
trol capability, it is recommended that the 97.7 percentile estimates of solar and gtxj 
magnetic activity be used. This then permits the design of a statistically conservative 
spacecraft lifetime at a given orbital altitude(s). It can Ix' interpi*etcd to mean the 
spacecraft will have a 97.7 percent chance of (1) remaining in orbit for at least the 
number of months calculateil by the orbital lifetime scheme and/or (2) maintaining the 
specified ix>ntrol capability. The determination sl,ould be made based upon the most 
current long range statistical estimate of the solar flux and geomagnetic index consist 
ent with the critical pi'oject development decision time points prior to launch of the 
spacecraft . 

For spacecraft already in orbit, or for some prolaunch orbi'.al altitude tradeoff 
analyses, the 97.7. 50, and 2.3 percentile values may be used to conq>ute the respec 
tive expected future lifetimes and control capability and then a stati.^'tical interpola 
tion made as a first ostima'e of tlx' probability that (1) the spacecraft will remain in 
orbit, or above a given orbital altitude, for any selected time period or (2) the 
six'cified control capability will l)e maintained for the required lime in orbit. 

After lauiu'h. the statements on remaining lifetime should always bo given in 
terms of the probability that a spacecraft's lifeiinx) will oiiual or exceed a given dale 
or dates depending on program managemeiU requirements. 

Since orbital lifotinx' programs and wntrol analyses require a specific date to 
associate with the solar flux F^^ ^ and geomagnetic A^^ cstMimtions to compute a 

corresponding atmospheric density, the data |X)ints should bo identified with the 
midway point of the given month. 
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It should be recognized and remembered that all statements of spacecraft orbital 
lifetime and control capabilities made using inputs of expected solar activity to the 
upper atmosphere models, must be interpreted in terms of statistical probabilities. 
Figure 46 provides a description of the various inputs necessary to compute a space* 
craft’s orbital lifetime. 



Figure 46. Solar predictions and spacecraft orbital lifetime. 


IX. CONCLUSIONS 


The results of this study strongly suggest that better predictions are possible, 
in a chi-square sense, by "lining up” the maximums (or minimums, or both). 

Evidence has also been presented which supports the existence of an aperiodic 
variation in the periods as well as in the amplitudes. 


34 





APPENDIX A 


ORIGINAL PAGE IS 
OF POOR QUALITY 


Least Squares Development 


The expected value of R (solar flux) at k points beyond m, where m is the 
number of points known (or used as known) in the current cycle and k is the 
number of points beyond m, is: 


m 


^ + 2 (C . AR . + b , Aft .) , 

m+k,j m+k q.k q,j q,k q,j » 


(A-1) 


where 


N 


*m+k N *^m+k,j ; 

^^m+k,j ~ ^m+k,j ^m+k 


N = Number of known past cycles 


(A-2) 

(A-3) 


Then Eq. (A-1) may be written as 
m 

AR = I (C . AR . + b .Aft .), 

m+k.J q,k q,j q.k q.j' 

Then 


(A- 4) 


Aft . = oTV- [AR ..1 • - AR , .] 
q,] 2At^. q+l,j q-l,J 

m 

AR = E (C , AR . + b , AR .) 

m+k,J q.k q.j q.k q.]^ 

The sum of the deviations squared is then 
N 

S = E 


3S 


3C. 


3s 

3b_, 


m 



AR - E 

m+k.] 


. AR . + b , 

,k q.j q.k 

N r 

m 


E AR - 

j=l 

E 

q=l 

C , AR . + b 
q.k q,] ( 

N r 

m 


jfi [ “‘“--■‘.i ■ 

E 

q=l 

C . AR . + b 

q.j q.l q 


AR 


(A -5) 


(A-6) 


(A -7) 


P.j (A-8) 


(A- 9) 
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This equated to zero gives the normal equations 
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N 


N m 


N m 


Note: 


^ AR , . AR^ . = Z Z C . AR , AR . + Z Z b , AR .AR . 
m+k.] p.] q,k q,j p.j q.k “"'q.] PO 

KT V, (A- 10) 

N N m N m 

^ ^^m+k 1 i = " ^ C n V ^ ^ Z b . AR , AR . . 

j=l P'3 i=i q=i q-J P>] j=i q=i Q.k q,] p ,j 

(A- 11) 


q,] 2At. q+l,j q-l,J 

P. . M 

At. = -2- ; P = s P- • 

2 P N j=l 3 


(A- 12) 
(A- 13) 


For one -point prediction, p = 1, m = 1. 


N 


N 


N 


Z AR , . AR . = Z C.. AR.. AR.. + Z b,, AR,. AR,. 
j=l k+l,j 1] ._j^ Ik 1] Ij Ik 1] 1] 


(A- 14) 


N 


N 


jfj ''Vl.i ''«p,j = jfj 4Ry ARy . ARy ARy . (A-15) 

Since the C's and b's are not summed, they may be taken outside the summation; 
i.e. , 


Cik ZARj. AR^. . bjj^ ZAR^. AR^ = IIAR^^^ . AR^ 


(A-16) 


^Ik ■*■ ^ik ^^^ik ^*^lj ” ^^^k+l,j ^^Ij 


(A-17) 


This can be written in vector and matrix notation if we let 



(A -18) 


(A-19) 
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(A- 20) 


-1 


Then AC = B -»• C = A B. 

For two -point prediction 

^^2k^®2]‘^*^lj (A-21) 

^^®2+k,j^*^2j ^ ^^2k^^2j^®2] ^^lk^®lj^®2j ^^2k^^2j^^2j (A-22) 

^^®2+k,j^^lj ^ * ^^2k^^2j^^lj ■*■ * ^^2k^^2j^^lj (A-23) 

^^*^2+k,j^^2j " ^^lk^^lj^^2j ^^2k^^2;^^2j (A-24) 

This may be written in vector and matrix notation if we let 



N 

, where I = Z 

j=l 


(A- 25) 



/zARjj 

ARy 

lAR^ 

ARy 

lAfty 

EAfty 

ARy 

A= 

E4Rl, 

*®2j 

EARj^. 

"®2j 

ZA^-. AR« 
1] 4 

ZAR2j 

''"21 


1 “«H 

Afty 

ZARo; 

Afty 

ZAfty ARj^j 

EASy 

Afty 


V““ii 


EiRy 

''*21 

zaR<. ARa; 

Ij 2j 

EAlly 

''*2i 



Then AC «= B -»■ C = A"^ B. 


(A-26) 


(A-27) 


(A-28) 
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APPENDIX B 

Cubic Connection From Predicted To Mean Cycle 


1) Select the last inflection point of the predicted cycle (tj, Fj) near the 

24-month point before the maximum. This should be done with at least a five-point 
numerical scheme for determining where the second derivative goes to zero. One 
such scheme is Sterling's [16] 



(B-1) 


2) Select the maximum (or minimum) on the mean cycle (t . F ) . This is 
normally input, but it can also be determined numerically. ™ 



(B-2) 


The second derivative should be zero at (tj, Fj), and the first derivative is 

zero at (t^, Fj„)* A cubic curve is the lowest degi*ee polynomial which can be 

determined that will go through the two points and have tliese properties at the two 
points . 

' The coefficients of the polynomial 

F = at^ + bt^ + ct + d (B-3) 


are to be determined from the following foiu* linear equations and will have the 
required four properties: 


6rtj + 2b = 0 


3at + 2bt + c = 0 
m m 


atj^ + btj + ctj + d = Fj 


+ bt„^ + ct„ + d = F„ 
m m m m 


I second derivative is zero at 
\ inflection point 

I first derivative is zero at 
I maximum or minimum 

I curve goes through the inflection 
I point (tj, Fj) 

I curve goes through the maximum 
\ (or minimum) point (t^^^, F^). 


(B-4) 

(B-5) 

(B-6) 

(B-7) 
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These four linear equations in the four coefficients are to be solved for a, b, c, and 
d. The solutions are: 



(B-8) 

(B-9) 

(B-10) 

(B-11) 


page fS 
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APPENDIX C 

Computer Program Listings 

INTC6CR TITLEU2) «PCYC ttPCV C2tCASE 

OmCNSION EPOCH <211 •PCRIOOI 211 tRt 132*2 11 tRNEAN<26*l t 
% RPRE0<26«1»YEAR<26A1 *R0eSV<26Al *cRRUP<26«l 1 *ERR0N 1 260 1 * 

% CH1PRS<2641*CH1TAB<26R1*06FREE<26R1 

INCLUDE PL0T1*LIST 
PtOTl pRoC 

nlHcNSION lTlT|.l< 12 1 « I TITL2 1 12 1 *I TITL3 < 12 1 * I TITL4| 121 

DIMENSION ICAPlll31»lCAP2<171*lCAP«m01*INsTR<221 

DATA INSTR/SoHFICHE AND ONE HARDCOPY. THANKS *1TP6H / 

DATA ITITL1/6HYEAR *11*6H / 

data ITITL2/18HS0LAR FLUX F19.7 *9*6H / 

DATA ItItL3/18HDECREES OF FREEDOM, 9*6H / 

data ITiTL 4/24HCUMULATIVC CHl SQUARE *8*6H / 

DATA ICAP1/36H*«* deROTCS ACTUAL *•* DENOTES NOMIN* 

$ «2KAI •♦•*•-• DENOTE *2 SIGMAS. -2 SIGpAS / 

data ICAP4/3CH*** DENOTES THEORETICAL DE, 

% 30HN0TES EXPERIMENTAL (PEARS0N*S1/ 

DATA IASTR/40/tI00T/35/»IPLUS/3<l/ .IMlNUS/sj/ 

1500 FoRMATIIS,* COEFFICIENT LINEAR MODEL USING CYCLES •,I?» 

% * THROUGH *»I2** PREDICTING CYCLES **I2»* AND **I2, 

% * FROM POINT **13*4X1 

CALL lOENTl 10$*INSTR1 

END 

READ! 5*50001 TITLE 
►JEAC«5*5C051 H,N 
EAOtStSClDl CEPOCHIJl ,J=1*N1 
70 5 J=1*N 

ITEaDIS.SOISI PERlOOtJl 
REA0(5*SC201 <R<I*J1*1=1*M1 
5 CONTINUE 

REA0(5*S025I NCASES 

DO 9C0 CASErl.NCXSES 

REA0<S*5C301 NUl* NU2,MU*MU2 .ALPHA 

PCYCl : HiiZ I 

PCYC2 = nU2 ♦ 2 

DO 10 K=i*M 

ROBSVIKl r R(K*FCYC11 

ROBSVIH^Nl = R<K*PCYC?1 

YEARIKI : EPOCHIPCYCll ♦ <K -1 1 ♦<PERIOD IPcYc 1 1 /Ml 
YCAR<M*K| : EP0CHIPCYC21 ♦ (K- 11* <PEPIOO< PC YC? 1/ Ml 
10 CONTINUE 

CALL LiNCMC<R*n*N.NU1.NU2.Mu*Mu 2*RNCAN,RPRE0»CRRUP.ERR0N1 
INCLUDE CULIST 
Cl PROC 

CALL PCHHAX<RPRE0*CRRUP*ERR0N*M,INFL«NPTS1 
MM : M ♦ NPTS 

END 

l*RITE( 6,65001 TITLE ,NU,NUl, NU2 .MUT.PCYCj *PC YC2 
«RlTE<6,6n35l ALPHA 
OO 20 K=1*MU2 

WRITE <6 *60101 YEARtKI *RObSV ( k> *RmEAN<K 1 
20 CONTINUE 

MU2P1 : HU2 ♦ I 
CHIRR = 0. 
include C2*LIST 
C2 PrOC 
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00 70 K:»».U2Pl,f<r' 

IF(K .67. INFL) GO TO 3C 
oev = RPRtO(K) • RHFAN(K) 

33 ChiPR = CHIpR ♦ (rPrEOIKI - ROBSV IK JI**2/RPPE0 «K > 

CHIPRSJK) = CHIPR 
NOF = K - MUZ 
OGFRCEIK) = KOF 

CHITAB(K) = CHIN(ALPHA,N0F,S40) 

60 TO SC 

40 CHITABIK) = 0.5*l SORT! 2,*D6FREEIK I - 1.1 ♦ l.Q602)*«2 
50 IFIK .GT. INFU 60 TO 60 

4RlTE(6t6:i5> VEARIK) ,ERRUP(K ) ,ROBS V IK ) ,RPR EC ( K ) ,ERRDH (K ) * 
« RMEA|^(K),OEytCHIPRS(K)»CHITASIK)»NOF 
CO TO 70 

6q write ( 6t6C20) V EAR (K 1 , ERRUP (K ) .ROBSV (k ) f RPREO I K ) »E rRONIK ) « 
%CHIPR5(K ItCHlTAC(K) ,NOF 
IF|K .LT. H) WRITEt6»6C25) 

70 CONTINUE 
END 

INCLUDE PL0T2,LIST 
PL0T2 PROC 

XMIN r FLOATilFIXIVEARllI - 1.)) 

XMAX : FlCATIIFIXIYEARI2*M) ♦ I,)) 

XMID : IXMIN * XMAX)/2, 

NP = MM - MU2PI ♦ 1 
CALL SETMIV|24, 24«120t24) 

CALL GRIOIV (3.xMIN,XMAXfO., 300.,. 5, 5.t.2»*2 *-4t-4»4t3) 
call xsclviixmik.ixmin.ierr) 

C«LL YSCLVK ISO., lYpIO.IERR I 
IX : IXMIN - 4J 
ly = IyMIO ♦ 9j 

call APRNtW»u.-12,16,ITITl2,IX ,IY ) 

CALL XSCLV1(XM1C,IXhI0,IErR) 

CALL YSCLVl(C.,IYMIN,IcRR) 

IX : IXMID - 12 
lY : lYMlN - 24 
CALL PRlNTV<4tITITLl,lX,IY) 

IX = IXMIO - 312 
lY : lY - 25 

CALL PRINTV<78» ICAPl.IX.IYI 
IX : IXMIO - 288 
lY ; lY - 15 

CALL PRINTV(72, title, lx, lYI 

EnCcOE ( 1C2,1C00 tlCAP2) Mu ,NUl , NU2 ,PCYC I ,PCY C 2 ,MU2 
IX : iXMjO - 408 
lY = lY - 15 

CALL PRINTV(102«ICAP2,1X,IY) 

CALL CHARSZtlS) 

CALL PL0TLN|-2*M,YEAR,R0BSV,IASTR ) 

call PL0TLN<-NP,YEARjMU2PU ,RPRE0<MU2P1I,ID0T» 

call PL0TlN('NP,YEARIMu2P1) , ERRUP (MU2P1) ,IPLUS) 

Call PLCTLNI -NP ,YEAP(MU2p1 > ,ERRDN(MU2P1I,IMINUS) 

CALL CHARSZaO) 

VMAX = AMAXO(CHIPRS(MM),CHITab(HN)I ♦ 10. 

YMID = YKAX/2. 

call DxDYVI2,0. ,YMAX,DY,MY, JY,NY,8. ,IERR) 

call GrIOIVI 3,0.,254.,U.,YMAX,5.,OY,-4,-MY, -4,-JY,3,NYJ 
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Call XSCLVKC.tlXMlNtlLf^R) 

CALL YSClVK VHID.IVMlotlERR I 
IX r IXMIN - «SY*2I*8 
lY = lYHlO ♦ 13s 

CALL APRNTV(0t*12«21s ITITL4 tlX tlY ) 

CALL XSCLVlU32.f IXMIO.ICRR I 
CALL YSCLV1(U«, lYMINtltRR) 

IX = IXMJ'> - 72 
lY : lYKIN - 24 
call PRIi^TVaS, 1T1 tL3,1X,Iy » 

IX = IXPIC - 2hC 
lY r IV - 25 

CALL PRlNTV(60tICAP<l,IX»IYI 
IX = IXMIO - 28« 

IV r IV - 15 

call PRlNTVI72,TITLe,lX,lY» 

IX : IXHIO - 43€ 

IV : IV - 15 

CALL PRlNTV(102*ICAP2*iXtIY ) 
call CHARS2(18I 

CALL PLOTLNC'NP •OGFRCE (HU 2 P 1 1 1 CHI TAB< MU2P1 1 » IaSTRI 
CALL PLOILNI -NP ,06FREC (HU2P I ) , CHI PR5 ( MU2Pi ) « » 
call CHARSZU4I 

END 

9oO CONTINUE 

CALL ENDJOB 
STOP 

SCOO F0RNAT(t>A6l 
5005 F0RHAT(2I4) 

5010 F0RHAT(8F9.3I 
5015 F0RhAT(4X,F5.1I 
5C20 F0RNAT(I2F6.n 
5C25 F0RMAT(l4l 
5030 F0RNAT(4I4,F4.2 I 

6000 FORHATt 1H1,12A6/1HO»I3,1X,*COEFFICIENT LINEAR HOOEL USING CYCLES*. 
% 1X,I2.1x»*THR0UGH*,1X.I2, 1X,»FR0H PO INT* . I X. 13 .IX .* PREDICTING *, 

S *CYClCS*»1X*I2.IX.*AnO*.1x*12) 

6C05 F0fiHATUHO,Tl04 .‘PEARSOn* *S * , 5X . * CHl-SOUARE * /6 x. *0ATE* , BX , * ^2 SIGH 
$AS* ,5X,*OBSERW£0*,6X,*PREoIcTED*. 5X,*-2 SIGH AS •,4X, ‘MEAN CyCLE*. 
«5X, •DEVIATION*, 5X,* CHI-SQUARE* »4X,*SIGNIFlcA NCE’ /IX# T102. 
S*(CUMULATIVE|#,4X,*LEVEL=SF4.2//I 
6C1 q F0RHAT<2X,Fe.3, 24X.F5.1.37X ,F5.1l 

6015 fORMATi2X.F8.3, lOX . 51 FS. l,9x> t F5. 1# 2X* 2< SX. F 1C.3 1 , • f • . 13. * I * 1 
^020 F0RHAT(2X,F8.3» 13X , 4( F5. 1 ,9X ) ,21X ,2 ( 5X ,F10. 3 I , * ( * . 1 3. * ) * I 
6025 FORHATliH-#,***! 

END 


ECAST-A 
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SUBROUTINE lInC kC f R »H «N»Nu 1 *NU2»NU t MU 2 •RHC;tN iRPTEU, CRRUR« ERRONl 
INTEGER P»Q 

DIMENSION RlMtN)»RMEAN(n,RPREo<l itERRUPllI ttRRDN< II* 

S OB(132,21),AI lO.lOl* AiNVf IQ » ID I » APRINE ( IQQ I ,a U32 t iO I r 
S Cll32»IQ)t JDUH|iai,0ETAi2) 
equivalence (A( l,l),AINV(l,ll,APRlMEtlll 
NUl = HU2 * NU •» 1 
NU : NU2 - NUl * 1 
MUPRIM = M - MU2 
DO 15 I=1*M 
RHEANtl) = Q. 

DO 10 J=NU1»NU2 

RmEAnII) = RMEAN(I) ♦ R(I,J) 

10 CONTINUE 

RMEANII) = RMEaN(I)/NU 
15 CONTINUE 
DO 25 1=1, M 
DO 2C J=NU1,N 

OR(I,J) = RiI,J| - RMEANII) 

20 continue 

25 continue 

DO RO K=1,MU 
00 35 0=1, MU 
A|K,Q) : 0. 
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CO 1C J=NUJlyNU2 

AIK^Q) = •(KfO) « OR|Hi;l4K*l, JI«OR(HUl4Q-l» J) 

3P CONTlN'Jf 
35 CONTINUE 
4 q coNrxNue 

00 55 P=l,r»UPRIM 
CU SO 0-1, HU 
BCP.UJ : 0. 

Do 45 J=KUl,Nu2 

BCP.O) = BIP,0» ♦ PR(HU24p, Mul»0-l,J> 
i«5 COKITINUC 
50 CONTINUE 
55 CONTINUE 
H = 0 

CO 65 U=i»HU 
DO 6C I=IfHU 
K r K ♦ I 
APRIMf (K > ; A(I ,J) 

60 COi^VZnUC 
65 continue 

CeTAm = 3. 

0ETAI2) : Q» 

CALL GJR(A,HUt»U»HU,HU»S999,J0UK,0ETA) 

IFCDETAIli ,EQ. 0.) GO TO 999 
K = NU«HU * 1 
qO 75 J — HSJflf*! 

00 ?0 Uftb,!,-! 

K = K « j 

AiNVSTtJ) = APRIMEIKI 
70 CONTINUE 
75 CONTINUE 

DO 90 r = l,MUPRI?' 

00 85 Krl^HU 
C(P»K) r 0. 

PO 0=it«u 

C<p,K) -■ C»P,K) 4 8«PtO)*AINVCO,K) 

80 continue 

85 CONTINUE 
90 CONTINUE 

DO IqS P=1,h 

IF(P .6T. HU2> GO TO 9‘ 

RPREOIP) = RIP.NUZ^l) 

60 T 3 ijjS 

95 RPREd»P) = 0. 

DO ItO K-1,MU 

RPRcD<P» : RPRCCIP) ♦ Cl P-HU2 , K » *DR < MUl 4 k-1 , nU24 1 ) 
100 CONTINUE 

RPREOlPl = RPRtOlP) ♦ RMEANIP) 

105 CONTINUE 

MU2PI = MU2 ♦ 1 
Do J20 P=Mu2Pl,P 
STGnA : q. 

00 115 j-NU1»nD2 

SIGMA : sigma ^ (R(P,JJ - RPREDf p ) »*«'2 
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11$ CONTINUE 

SIGMA = SQRT«SI6MA/INU-MU^in 
CRRUPIP) = RPREDIPI ♦ 2.*SIGHA 
ERR0N(Pl = RPREDIPI • 2»*SICMA 
IFIERRONlPl .LT , O#! ERRON«PI 
120 CONTINUE 
RETURN 

999 STOP SINGLR 

END 
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